Detecting Multimodal AUC distributions in HNSCC¶
We currently assign sensitivities by tail-assignment, however, there may be cases with distinctive multimodal distributions that are indicative of different response models.
We'll begin by testing the method on beatAML data where we have much greater number of observations, and therefore more well defined distributions.
$$ p(x) = \frac{1}{\sqrt{ 2 \pi \sigma^2 }} e^{ - \frac{ (x - \mu)^2 } {2 \sigma^2} } $$
In [1]:
%%javascript
IPython.OutputArea.prototype._should_scroll = function(lines) {
return false;
}
In [2]:
import pandas as pd
import numpy as np
from sklearn.mixture import GMM
import seaborn as sbn
from matplotlib import pyplot as plt
import warnings
# ------ for mut models ---------
from sklearn.linear_model import ElasticNet
from sklearn.linear_model import LogisticRegression
from sklearn.manifold import TSNE
import statsmodels.api as sm
from sklearn import metrics
warnings.filterwarnings("ignore")
In [12]:
# import jtplot module in notebook
from jupyterthemes import jtplot
# choose which theme to inherit plotting style from
# onedork | grade3 | oceans16 | chesterish | monokai | solarizedl | solarizedd
jtplot.style(theme='onedork')
In [4]:
aml_aucs = pd.read_csv('./../beatAML_data/beatAML_aucs.csv')[['inhibitor','lab_id','auc']].drop_duplicates()
#aml_aucs = pd.read_csv('./../data/beatAML_AUCs_subset.csv')
In [5]:
aml_aucs.head()
Out[5]:
In [6]:
def add_normal_plot(mu, s, weight, c, ax, auc_max=300):
'''
mu = mean
s = standard deviation
ax = matplotlib axes to add to
c = color {[r b g c y ...]} <str>
'''
x = np.arange(0, auc_max, 1)
y = (1/(2*np.pi*s**2)**(0.5))*np.exp( - (x-mu)**2/(2*s**2) ) * weight
ax = ax.plot(x,y, color=c, label='mean: %.1f, std: %.1f' %(mu, s))
def get_color():
'''
'''
for c in ['r','b','g','c','y']:
yield c
In [7]:
gen = get_color()
gen.__next__()
Out[7]:
In [8]:
def test_multimodal_fits(X, ntests=10, kmax=5, plot=True, inhib=None, override_k=False):
'''
'''
res = {x:[] for x in ['k', 'aic', 'bic']}
for k in range(1,kmax):
for i in range(ntests):
#print('k: %d' %k)
gmm = GMM(n_components=k, n_init=1)
gmm.fit(X)
res['k'].append( k )
res['aic'].append( gmm.aic(X) )
res['bic'].append( gmm.bic(X) )
res = pd.DataFrame( res )
best_k = res[res.bic == np.min(res.bic)].k.unique()[0] if not override_k else override_k
gmm_best = GMM(n_components=best_k,n_init=20)
gmm_best.fit(X)
P = gmm_best.predict(X)
if plot:
nbins = 50
bin_ = np.arange(0,np.max(X),np.max(X)/nbins)
scalar_to_make_pretty = 0.25 # since our fitted Gaussians are normalized to their weights, they appear smaller
weights_ = scalar_to_make_pretty/len(X)
f, axs = plt.subplots(1,3,figsize=(15,5))
sbn.distplot(X, bins=bin_, ax=axs[0]).set_title('AUC distribution')
sbn.scatterplot(x='k', y='bic', alpha=0.3, data=res, ax=axs[1]).set_title('BIC vs K')
clas = 0
for weight, mean, covars, c in zip(gmm_best.weights_, gmm_best.means_, gmm_best.covars_, get_color()):
sbn.distplot(X[P==clas], bins=bin_, kde=False, color=c, ax=axs[2], label='AUC', hist_kws={'weights': [weights_]*len(X[P==clas])})
add_normal_plot(mean[0], (covars[0])**0.5, weight, c, axs[2], auc_max=np.max(X))
clas+=1
axs[2].set_title('Optimal GMM fit')
plt.legend()
plt.suptitle(inhib)
print('Number of assays (aucs): %d' %len(X))
print('Optimal K: %d [BIC=%.1f]' %(best_k, np.min(res.bic)))
print('GMM fit:\n\tMixture Weights: %r\n\tMeans: %r\n\tVariances: %r' %(gmm_best.weights_.ravel(), gmm_best.means_.ravel(), gmm_best.covars_.ravel()))
print('Class counts: %r' %['class %d: %d' %(cl, len(X[P==cl])) for cl in list(set(P))])
plt.show()
return P
In [13]:
for inhib in aml_aucs.inhibitor.unique():
inhib_dat = aml_aucs[aml_aucs.inhibitor == inhib]
if inhib_dat.shape[0] > 450:
print('------------------------------------------------')
print('Inhibitor: %s' %inhib)
print('------------------------------------------------')
AUCS = inhib_dat.auc.values.reshape(-1,1)
test_multimodal_fits(AUCS, ntests=10, kmax=6, inhib=inhib, plot=True)
Save subset of okay'd data¶
In [10]:
okay_data = ['YM-155', 'Vandetanib (ZD6474)', 'Trametinib (GSK1120212)', 'Sunitinib','Sorafenib', 'Selumetinib (AZD6244)', 'Nilotinib', 'JAK Inhibitor I', 'Pazopanib (GW786034)', 'Elesclomol', 'Dasatinib']
dat = aml_aucs[aml_aucs.inhibitor.isin(okay_data)]
dat.to_csv('./../data/beatAML_AUCs_subset.csv')
Test HNSCC distribution¶
The next step is to ask, can we model these drug responses across cancer types? As in, are the distributions that are well defined by the beatAML similarly define the HNSCC distribution?
We will try to test this by setting our null hypothesis as:
Let X be the underlying set distribution such that: $$ AUC_{beatAML} \in X $$
$$ H_0: AUC_{HNSCC} \in X $$ and $$ H_1: AUC_{HNSCC} \notin X $$
We will use permutation testing, and reject cases where alpha=0.05.
In [14]:
HNSCC_all = pd.read_csv('./../beatAML_data/HNSCC_all_functional_data.csv')
HNSCC_auc = HNSCC_all[['lab_id','inhibitor','auc','call']].drop_duplicates()
HNSCC_auc.head()
Out[14]:
In [15]:
shared_drugs = list( set(HNSCC_all.inhibitor) & set(aml_aucs.inhibitor) )
shared_drugs
Out[15]:
In [19]:
def permutation_test(x, y, n=1e5, verbose=True, return_prob=True, alpha=0.05):
'''
two tailed permutation test to see if y is apart of same model as x
returns the probability that y is not drawn from the same underlying model (or True if prob > 0.05 and return_prob=True)
if pval > alpha -> x,y may not have been drawn from same distribution (alpha probability not drawn)
H0: Distributions are the same
H1: the distributions are the not the same
but we want 95% confidence that the distributions are the same... so,
'''
x = x.ravel()
y = y.ravel()
z = np.append(x,y)
T = np.abs( np.mean(x) - np.mean(y) )
permutation_means = []
for i in range(int(n)):
if (i % 333 == 0) and verbose: print('Running permutations...[%.2f%%]' %(i/n*100), end='\r')
perms = np.random.permutation(z)
perm_y = perms[:len(y)]
perm_x = perms[len(y):]
permutation_means.append( np.abs(np.mean(perm_y) - np.mean(perm_x)) )
print()
print(np.sum(1* np.array(permutation_means) > T))
pval = np.sum( 1*np.array(permutation_means) > T ) / len( permutation_means )
return pval if return_prob else pval >= alpha
def test_distribution_similarity(aml, hnscc, inhib=None, verbose=True, perms=1e5):
'''
'''
if verbose: print('---------------------------------------------------------')
if verbose: print('XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX')
if verbose: print(inhib)
if verbose: print('XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX')
if verbose: print('---------------------------------------------------------')
if verbose: print()
print('HNSCC assays: %d' %hnscc.shape[0])
print('AML assays: %d' %aml.shape[0])
test_multimodal_fits(aml, ntests=10, kmax=6, inhib=inhib, plot=True)
test_multimodal_fits(hnscc, ntests=10, kmax=6, inhib=inhib, plot=True)
test_multimodal_fits(np.append(aml.ravel(),hnscc.ravel()).reshape(-1,1), ntests=10, kmax=6, inhib=inhib + ' hnscc-aml', plot=True)
nbins=50
bin_ = np.arange(0,300,300/nbins)
plt.figure()
plt.hist(aml_,color='blue',label='AML_AUC',normed=True,alpha=0.4,bins=bin_)
plt.hist(hnscc_,color='red',label='HNSCC_AUC',normed=True,alpha=0.4,bins=bin_)
plt.legend()
plt.show()
pval = permutation_test(aml_, hnscc_, n=perms, verbose=True, return_prob=True)
if verbose: print('Dasatinib permutation test p-value: %f' %pval)
In [21]:
for inhib in shared_drugs:
aml_ = aml_aucs[aml_aucs.inhibitor == inhib].drop_duplicates().auc.values.reshape(-1,1)
hnscc_ = 100*HNSCC_all[HNSCC_all.inhibitor == inhib].auc.drop_duplicates().values.reshape(-1,1) #hnscc is weighted by 1/100
if aml_.shape[0] + hnscc_.shape[0] > 300:
test_distribution_similarity(aml_, hnscc_, inhib=inhib, perms=1e5)
In [22]:
AML_Dasatinib = aml_aucs[aml_aucs['inhibitor'] == 'Sorafenib']
print('Number of AML Sorafenib assays: %d' %len(AML_Dasatinib))
AML_Dasatinib.head()
Out[22]:
In [23]:
HNSCC_Dasatinib = HNSCC_auc[HNSCC_auc['inhibitor'] == 'Sorafenib'].drop_duplicates()
print('Number of HNSCC Sorafenib assays: %d' %len(HNSCC_Dasatinib))
HNSCC_Dasatinib.head()
Out[23]:
In [24]:
print('------------------------------------------------')
print('Inhibitor: %s' %'AML_Dasatinib')
print('------------------------------------------------')
aml_= AML_Dasatinib.auc.values.reshape(-1,1)
c = test_multimodal_fits(aml_, ntests=10, kmax=6, inhib=inhib, plot=True)
print('------------------------------------------------')
print('Inhibitor: %s' %'HNSCC_Dasatinib')
print('------------------------------------------------')
hnscc_=100*HNSCC_Dasatinib.auc.values.reshape(-1,1)
c = test_multimodal_fits(hnscc_, ntests=10, kmax=6, inhib=inhib, plot=True)
In [25]:
nbins=50
bin_ = np.arange(0,300,300/nbins)
plt.figure()
plt.hist(aml_,color='blue',label='AML_AUC',normed=True,alpha=0.4,bins=bin_)
plt.hist(hnscc_,color='red',label='HNSCC_AUC',normed=True,alpha=0.4,bins=bin_)
plt.legend()
plt.show()
In [26]:
pval = permutation_test(aml_, hnscc_, n=1e6, verbose=True, return_prob=True)
print('Dasatinib permutation test p-value: %f' %pval)
Drug response associated mutations¶
In [27]:
aml_var_all = pd.read_csv('./../data/vizome_beatAML_variants.csv', skiprows=[0,1,2])
#[aml_var_all.tumor_only == True]
aml_var = aml_var_all.dropna()[['sample_id', 'gene', 'chr','start','end','ref','alt','variant','sift']]
aml_var = aml_var.rename(columns={'sample_id':'lab_id'})
print(aml_var.head())
print(aml_var[['chr','start','end','ref','alt']].head())
print('Number of genes mutated: %d' %len(aml_var.gene.unique()))
print('Number of unique mutations: %d' %aml_var[['chr','start','end','ref','alt']].drop_duplicates().shape[0])
mut_cnt = aml_var.groupby('gene').count()['lab_id'].sort_values()
mult_patients = mut_cnt.index[mut_cnt > 20]
#aml_variants = aml_var[aml_var['lab_id'].isin(mult_patients)]
shared_mutations = [x for x in mult_patients]
print(shared_mutations)
print(aml_var.groupby('gene').count()['lab_id'].sort_values())
aml_var = aml_var[aml_var['gene'].isin(shared_mutations)]
aml_var.head()
In [ ]:
# convert variants to matrix format
genes = pd.Series( aml_var.gene.unique() )
res = []
for lab_id in aml_var.lab_id.unique():
D = aml_var[aml_var['lab_id'] == lab_id]
x = genes.isin( D['gene'] ).tolist()
res.append( [lab_id] + [str(int(xx)) for xx in x] )
res = pd.DataFrame( res, columns=['lab_id'] + genes.tolist() )
res.head()
MAT_VAR = res
print('matrix variant shape: %s' %str(res.shape))
print(MAT_VAR.head())
In [ ]:
print( aml_aucs.head() )
aml_aucs['lab_id']=aml_aucs['lab_id'].astype(str)
MAT_VAR['lab_id']=MAT_VAR['lab_id'].astype(str)
aml_auc_w_var = aml_aucs.merge(right=MAT_VAR, how='inner', on=['lab_id'])
aml_auc_w_var.head(5)
In [ ]:
Dasatinib = aml_auc_w_var[aml_auc_w_var['inhibitor'] == 'Dasatinib']
print('Dasatinib auc and variants shape: %s' %str(Dasatinib.shape))
aucs = Dasatinib.auc.values.reshape(-1,1)
X = Dasatinib[genes]
Y = test_multimodal_fits(aucs, inhib='Dasatinib', override_k=2)
X_embedded = TSNE(n_components=2).fit_transform(X)
df = pd.DataFrame( {'x1':X_embedded[:,0], 'x2':X_embedded[:,1], 'label':Y} )
plt.figure()
sbn.scatterplot(x='x1', y='x2', hue='label', data=df)
plt.show()
In [ ]:
reg = ElasticNet().fit(X, Y)
yhat = reg.predict(X)
pred = (yhat > 0.5)*1
acc = sum([a == b for a,b in zip(pred, Y)]*1) / len(Y)
print('Elastic Net accuracy: %.2f' %acc)
fpr, tpr, thresholds = metrics.roc_curve(Y, yhat, pos_label=1)
EN_auc = metrics.auc(fpr, tpr)
print('Elastic Net AUC: %.2f' %EN_auc)
plt.figure()
plt.plot(fpr,tpr,'r--')
plt.show()
logit = LogisticRegression().fit(X,Y)
yhat = logit.predict(X)
pred = (yhat > 0.5)*1
print(pred)
acc = sum([a == b for a,b in zip(pred, Y)]*1) / len(Y)
print('Logistic Regression accuracy: %.2f' %acc)
fpr, tpr, thresholds = metrics.roc_curve(Y, yhat, pos_label=1)
EN_auc = metrics.auc(fpr, tpr)
print('Logistic Regression AUC: %.2f' %EN_auc)
plt.figure()
plt.plot(fpr,tpr,'r--')
plt.show()
print('---------------------------------------------------')
print('STATSMODEL LOGIT')
print('---------------------------------------------------')
X2 = X.values.astype(int)
logit_sm = sm.Logit(Y, sm.add_constant(X2)).fit(method='lbfgs', maxiter=500)
yhat = logit_sm.predict(sm.add_constant(X2))
cutoff = 0.9
pred = (yhat > cutoff)*1
acc = sum([a == b for a,b in zip(pred, Y)]*1) / len(Y)
print('Logistic Regression accuracy: %.2f' %acc)
fpr, tpr, thresholds = metrics.roc_curve(Y, yhat, pos_label=1)
EN_auc = metrics.auc(fpr, tpr)
print('Logistic Regression AUC: %.2f' %EN_auc)
plt.figure()
plt.plot(fpr,tpr,'r--')
plt.show()
print( logit_sm.summary() )
#print(yhat)
#build_RF_variant_model(X,Y)